% MLE estimation of parameters
% modified for RAND round 1 revision on June 25, 2022
% RUN this on RIVANNA on Interactive mode. 
% Runs in within 2 hours. 

rng(5758934)

data0 = open('OilData.mat');
dataQ = table2array(data0.OilData.Production)*30/1000;     % dataQ is in million barrels per month
dataQ_names = data0.OilData.Production.Properties.VariableNames;

dataP = data0.OilData.Price;
% Deflator (2019 Q4 =100): EXPLAIN IT IN THE PAPER
T = size(dataQ,1); I = size(dataQ,2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% De-trending
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% detrending the data (x must be row vector)
tempfun = @(x,tdata) -exp(repmat(x(1:I),T,1)-x(2*I+1)*tdata)...
    + exp(repmat(x(I+1:2*I),T,1));
tdata = repmat((1:T)',1,I);

tempfun1 = @(x) sum(sum((dataQ - tempfun(x,tdata)).^2));
options = optimoptions('ga','MaxGenerations',250*(2*I+1));
tau_est = ga(tempfun1,2*I+1,[],[],[],[],zeros(1,2*I+1),[9^999*ones(1,2*I) 1],[],[],options);

%% IMPORTANT: hereafter we use detrended dataQ!
dataQ = dataQ + exp(repmat(tau_est(1:I),T,1)-tau_est(2*I+1)*tdata) ;

% compute means and std. deviations
dataQ_means = mean(dataQ);
dataQ_std = std(dataQ);
table0 = [dataQ_means' dataQ_std'];

%%% clustering in ng groups
% starting point of k-means algorithm (marked in graph)
%start_obs = [2 7 17 19]; ng = length(start_obs);

%[groupid, centroids] = kmeans(table0,ng,'MaxIter',10000,'Start',table0(start_obs,:));

% already did the Kmeans... then ised KMEANS + geograpphy to get the
% following

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% ESTIMATION %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 6 groups 
ng = 6; 

% new classification that combines K-means and Geography
groupid = [4;... %Angola
          1;... %UAE
          2;... %Canada
          2;... %China
          4;... %Algeria
          3;... %Ecudaor
          4;... %Egypt
          2;... %UK
          1;... %Iran
          1;... %Iraq
          1;... %Kuwait
          4;... %Libya
          3;... %Mexico
          4;... %Nigeria
          2;... %Norway
          1;... %Qatar
          5;... %Russia
          5;... %Saudi
          6;...  %USA
          3];  % Venezuela

data1 = [dataP dataQ];
conv_ma = (repmat(groupid,1,ng)==repmat(1:ng,I,1));

%we started with this
% param_ini = ...
%     [0.2/30    0.1/30   47.9592  513.2512...
%     0.3942    0.3851    0.2459  0.3851    0.2459  3.3491...
%     0.9347    0.6093    0.8824  0.6093    0.8824 10.9169...
%     0.0015    0.0039   18.8001];

%estimated several times and came up with the new initial guess
 
param_ini =[0.0215 ...
    0.0012 ...
   97.4799 ...
  551.1091 ...
    0.9410 ...
    1.0045 ...
    0.5492 ...
    0.8465 ...
    0.3648 ...
    0.1553 ...
    0.8727 ...
    0.8401 ...
    0.3372 ...
    0.5690 ...
    1.2447 ...
    0.9530 ...
    0.0001 ...
    0.0002 ...
   19.5021 ...
   57.1800];


tempres = estimate_fun(data1,conv_ma,I,ng,param_ini);
param_est = tempres(1:end-1);
llfvalue = tempres(end);


% beta_est = param_est_new(1);
% lambda_est = param_est_new(2);
% muU_est = param_est_new(3);
% sig2U_est = param_est_new(4);
% alv_est = param_est_new(5:ng_new+4);
% betv_est = param_est_new(ng_new+5:2*ng_new+4);
% a1_tilde_est = param_est_new(2*ng_new+5);
% a2_tilde_est = param_est_new(2*ng_new+6);
% barw_est = param_est_new(2*ng_new+7);
% ulow_est = param_est_new(end);

theta_mle      = [param_est(1), ... % beta
                  param_est(2), ... % lambda
                  param_est(3), ... % muU
                  param_est(4), ... % sig2U
                  param_est(5:ng+4),... %alv
                  param_est(ng+5:2*ng+4), ... % betv
                  param_est(2*ng+5), ... %a1-tilde
                  param_est(2*ng+6), ... %a2-tilde
                  param_est(2*ng+7), ... %bar-w
                  param_est(end), ... %ulow
                  ];


save('/home/ga6x/pd/cournot/RAND_revision_round1_ESTIMATES.mat',...
    'theta_mle','-v7.3')
